home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / blk_con.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  166 lines
  1. ;$Id: blk_con.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       BLK_CON
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes a "fast convolution" of a digital signal 
  11. ;       and an impulse-response sequence.
  12. ;
  13. ; CATEGORY:
  14. ;       Digital Signal Processing
  15. ;
  16. ; CALLING SEQUENCE:
  17. ;       Result = BLK_CON(Filter, Signal, B_length = B_length)
  18. ;
  19. ; INPUTS:
  20. ;       Filter = A P-element floating-point vector containing the impulse-
  21. ;                response sequence of the digital filter.
  22. ;       Signal = An N-element floating-point vector containing the discrete 
  23. ;                signal samples.
  24. ;                
  25. ; KEYWORD PARAMETERS:
  26. ;       B_length = (Block Length) An integer specifying the length of
  27. ;                  the subdivided signal segments. If this paramter is
  28. ;                  not specified, a near-optimal value is chosen by the 
  29. ;                  algorithm based upon the length of the impulse-response 
  30. ;                  sequence, P. If P is a value less than 11 or greater 
  31. ;                  than 377, then B_length must be specified.
  32. ;
  33. ; RESTRICTIONS:
  34. ;       1) The block length must be greater than the filter length.
  35. ;          B_length > P
  36. ;       2) The block length must be less than the number of signal samples.
  37. ;          B_length < N_elements(Signal)
  38. ;
  39. ; EXAMPLE:
  40. ;       Create an impulse-response sequence of length P = 32. 
  41. ;         filter = replicate(1.0, 32) ;Set all points to 1.0
  42. ;         filter(2*indgen(16)) = 0.5  ;Set even points to 0.5
  43. ;
  44. ;       Create a sampled signal with random noise.
  45. ;         signal = sin((findgen(1000)/35.0)^2.5)
  46. ;         noise  = (randomu(SEED,1000)-0.5)/2
  47. ;         signal = signal + noise
  48. ;
  49. ;       Convolve the filter and signal using block convolution.
  50. ;         result = BLK_CON(filter, signal)
  51. ;
  52. ; PROCEDURE:
  53. ;       Implementation of the "overlap-save" method in the frequency domain.
  54. ;       The discrete signal samples are divided into overlapping segments of
  55. ;       a length specified by the parameter B_length. B_length may be supplied
  56. ;       by the user as an optional keyword parameter or determined by the
  57. ;       algorithm to a near-optimal value. Each input segment consists of P-1
  58. ;       samples from the previous input segment and (B_length-P+1) new signal
  59. ;       samples, where P is the length of the filter. Each of these segments
  60. ;       is processed in the frequency domain and then 'reassembled' in the
  61. ;       time domain. The first and last input segments are handled differently.
  62. ;       The result is an N-element floating-point vector containing the 
  63. ;       filtered signal.
  64. ;
  65. ; REFERENCE:
  66. ;       Oppenheim, A.V. and Schafer, R.W.
  67. ;       DIGITAL SIGNAL PROCESSING
  68. ;       Prentice-Hall, 1975
  69. ;
  70. ; MODIFICATION HISTORY:
  71. ;           Written by:  GGS, RSI, May 1993
  72. ;           Modified:    GGS, RSI, June 1994
  73. ;                        Added long indexing into vectors. Minor changes in 
  74. ;                        the use of intermediate variables reduces memory
  75. ;                        allocation in certain instances. Made slight changes
  76. ;                        to the documentation header.
  77. ;-
  78.  
  79. function BLK_CON, Filter, Signal, B_length = B_length
  80.  
  81.   ;Return to caller if an error occurs.
  82.   on_error, 2 
  83.  
  84.   ;Block length is based upon filter length.
  85.   p = n_elements(filter)
  86.   if n_elements(b_length) eq 0 then begin
  87.     if (p lt 11) then stop, $
  88.       'Block length must be specified as a keyword parameter.' $
  89.     else if (p LE 17)  then blen = 64   $
  90.     else if (p LE 29)  then blen = 128  $
  91.     else if (p LE 52)  then blen = 256  $
  92.     else if (p LE 94)  then blen = 512  $
  93.     else if (p LE 171) then blen = 1024 $
  94.     else if (p LE 377) then blen = 2048 $
  95.     else stop, 'Block length must be specified as a keyword parameter.'
  96.   endif else blen = long(b_length) ;User specified block length.
  97.  
  98.   ;RESTRICTION:1
  99.   if p ge blen then stop, $
  100.     'Block length must be greater than the filter length.'
  101.  
  102.   ;Number of discrete signal samples.
  103.   ns = n_elements(signal) 
  104.  
  105.   ;RESTRICTION:2
  106.   if blen ge ns then stop, $
  107.     'Block length must be less than the number of signal samples.'
  108.  
  109.   ;Number of signal segments of length, blen-p+1
  110.   k = ns / (blen - P + 1L)
  111.  
  112.   ;Length of last signal segment.
  113.   rem = ns mod (blen - p + 1L)
  114.  
  115.   ;Represent the filter in frequency domain.
  116.   z_filter = fft([filter, fltarr(blen-p)], -1)
  117.  
  118.   ;The first segment is handled differently than the rest.
  119.   ;Forward zero-pad signal segment up to a length of blen.
  120.   lower = 0L
  121.   upper = blen - p
  122.   input = [fltarr(p-1), signal[lower:upper]]
  123.   
  124.   ;Allocate storage for the result.
  125.   result = fltarr(ns, /nozero)
  126.  
  127.   ;Complex point-wise multiplication and inverse transform.
  128.   output = float(fft((z_filter * fft(input, -1)), 1) *  blen)
  129.   ;The 'blen' term is a scale factor unique to IDL's FFT algorithm.
  130.   ;Discard the first p-1 points to obtain the first output segment.
  131.   result[0] = output[p-1:blen-1]
  132.  
  133.   ;Next place to store.
  134.   ip = blen - p + 1L
  135.  
  136.   ;Begin to process the subdivided signal segments.
  137.   loop  = 1L
  138.   lower = upper + 1L
  139.   upper = lower + blen - p  ;upper = lower+(blen-p+1)-1
  140.   while loop le k-1L do begin
  141.     input  = [input[blen-p+1:blen-1], signal[lower:upper]]
  142.     output = float(fft((z_filter * fft(input, -1) ),1) * blen)
  143.     result[ip] = output[p-1:blen-1] ;Next output segment.
  144.     ip = ip + blen - p + 1L
  145.     lower = upper + 1L
  146.     upper = lower + blen - p
  147.     loop  = loop + 1L
  148.   endwhile
  149.  
  150.   ;Last signal segment is of length, rem. If rem is not zero 
  151.   ;there is a signal segment that is handled differently.
  152.   if rem ne 0L then begin
  153.     output = float(fft((fft([filter, fltarr(p+rem-2)] ,-1) $
  154.              * fft([input[blen-p+1:blen-1], signal[lower:lower+rem-1], $
  155.                fltarr(p-1)], -1)), 1) * (2*p+rem-2))
  156.     ;(2*p+rem-2) is an FFT scale factor.
  157.     result[ip] = output[p-1:p+rem-2]
  158.   endif
  159.   return, result
  160.   ;Clean up intermediate variables.
  161.   z_filter = 0
  162.   input = 0
  163.   output = 0
  164.   result = 0
  165. end
  166.